%% Supplementary: TEM datasets.

%Load files. You may need to put the full directory name here. They are in    root\TEM
folderName = 'TEM\';
fileList = dir(strcat(folderName, '*.dm4'));
clear img;

[tags, img1] = dmread(strcat(folderName, fileList(1).name));
[tags, img2] = dmread(strcat(folderName, fileList(2).name));
[tags, img3] = dmread(strcat(folderName, fileList(3).name));
% analyze positions of atoms and see if they are consistent across layers. 

folderName = 'TEM\Simulations\';
fileList = dir(strcat(folderName, '*.dm4'));
clear img;

fileList([2:11]) = fileList([11, 2:10]);

[tags, simg1] = dmread(strcat(folderName, fileList(1).name));
[tags, simg2] = dmread(strcat(folderName, fileList(2).name));
[tags, simg3] = dmread(strcat(folderName, fileList(3).name));
[tags, simg4] = dmread(strcat(folderName, fileList(4).name));
[tags, simg5] = dmread(strcat(folderName, fileList(5).name));
[tags, simg6] = dmread(strcat(folderName, fileList(6).name));
[tags, simg7] = dmread(strcat(folderName, fileList(7).name));
[tags, simg8] = dmread(strcat(folderName, fileList(8).name));
[tags, simg9] = dmread(strcat(folderName, fileList(9).name));
[tags, simg10] = dmread(strcat(folderName, fileList(10).name));
[tags, simg11] = dmread(strcat(folderName, fileList(11).name));
scale = 166.775; %scale is 166.775 px / nm (measured from c direction).
x = [0:511] / scale;

%% Fig. S2D and Inset: Lineshape from Misorientation from Zone Axis. Fit widths for the inset. 
figure(1);
clf;
colors = jet(11);

cuts = nan(11, 512);
for i = 1:11
    str = strcat('data = simg', num2str(i), ';')
    eval(str); 
    cuts(i, :) = sum(data(225:235, :), 1);
    
    plot(x, cuts(i, :), 'color', colors(i, :)); hold all;

end

xlim([0.3, 0.8]);
xlabel('x (nm)');
ylabel('TEM line cut (arb)');
figure(2);
clf;
plot(x, cuts(1, :));
hold all;
plot(x, cuts(11, :));

% fit cuts.
lorentzFn = fittype('2*A/pi * (f./(4*(x-c).^2 + f.^2))');

for i = 1:11
    i
    %for each atomic site
    data = cuts(i, :) - mean(cuts(i, :));

    figure(20);
    clf;
    plot(x, data);
    hold all;
    
    %find the idx where the cuts cross the mean
    differences = sign(data(40:end-40));
    signs = differences(1:end-1) .* differences(2:end);
    idx = find(signs < 0) + 40;
    plot(x(idx), data(idx), 'o');
    
    %detect if the first cycle has a maximum in it. If not delete the first
    if data(round(mean([idx(1),idx(2)]))) < 0
        idx(1) = [];
    end
    
    for j = 1:floor(length(idx)/2)        
        elements = (idx(2*j-1)-10) : (idx(2*j)+10);
        midpt = mean(x(elements));
        y = data(elements)';
        y = y - min(y);
        [fitResult, fitgoodness] = fit(x(elements)', y, lorentzFn, 'StartPoint', [max(y), midpt, 0.03], 'Robust', 'off', 'Algorithm', 'Trust-Region', 'TolFun', 1e-10, 'TolX', 1e-10);

        atomsBot(i, j) = fitResult.c;
        atomsBotWidth(i, j) = fitResult.f;
         
        figure(55);
        clf;
        plot(x(elements), y, 'o');
        hold all;
        plot(x(elements), fitResult(x(elements)));
        xlim(fitResult.c + [-0.3, 0.3]);
        pause(0.05);
    end
    
    figure(20);
    plot(atomsBot(i, :), interp1(x, data, atomsBot(i, :)), 'x', 'color', [1,0,0]);
end

atomsBotWidth = mean(atomsBotWidth, 2);
% Fig. S2D Inset: Fitted widths of peaks.
figure(10);
clf;
plot((0:5:50) * 0.0572958, atomsBotWidth * 10, '^')
xlabel('Local Tilt Angle (degree)');
ylabel('FWHM (A)');
xlim([0, 1.5]);
ylim([0.08, 0.16] * 10);


%% Fit for positions of individual (columns of) atoms along BiO layers.
%take cuts along BiO layers.
list = [1944, 2030, 2118, 2206, 2292, 2380, 2468, 2555, 2643, 2730, 2817, 2904, 2992, 3080, 3167, 3255, 3342, 3430, 3517, 3605, 3692, 3780];
cutsBottom = zeros(length(list), 4096);
for i = 1:length(list)
    
    cutsBottom(i, :) = sum(img1(list(i)-1:list(i)+17, :), 1);
    imgX(list(i)-1:list(i)+17, :) = 0;
end

listB = [1934, 1848, 1760, 1673, 1585, 1499, 1410, 1323, 1235, 1148, 1061, 974, 884, 797, 710, 622, 535, 448, 360, 273, 186, 98];
cutsTop = zeros(length(listB), 4096);
for i = 1:length(listB)
    cutsTop(i, :) = sum(img1(listB(i)-14:listB(i)-5, :), 1);
    
    imgX(listB(i)-14:listB(i)-5, :) = 0;
end


%set up fitting function.
x = 74.31 * 0.9662319529 / 4096 * (0:4095);
lorentzFn = fittype('2*A/pi * (f./(4*(x-c).^2 + f.^2))');
atomsBot = NaN(22, 189);
atomsBotWidth = NaN(22, 189);
%for each BOTTOM layer,
for i = 1:22
    i
    data = cutsBottom(i, :) - mean(cutsBottom(i, :));

    figure(20);
    clf;
    plot(x, data);
    hold all;
    
    %find the idx where the cuts cross the mean
    differences = sign(data(40:end-40));
    signs = differences(1:end-1) .* differences(2:end);
    idx = find(signs < 0) + 40;
    plot(x(idx), data(idx), 'o');
    
    %detect if the first cycle has a maximum in it. If not delete the first
    if data(round(mean([idx(1),idx(2)]))) < 0
        idx(1) = [];
    end
    
    parfor j = 1:floor(length(idx)/2)        
        elements = (idx(2*j-1)-8) : (idx(2*j)+8);
        midpt = mean(x(elements));
        y = data(elements)';
        y = y - min(y);
        [fitResult, fitgoodness] = fit(x(elements)', y, lorentzFn, 'StartPoint', [max(y), midpt, 0.03], 'Robust', 'off', 'Algorithm', 'Trust-Region', 'TolFun', 1e-10, 'TolX', 1e-10);

        atomsBot(i, j) = fitResult.c;
        atomsBotWidth(i, j) = fitResult.f;
    end
    
    figure(20);
    plot(atomsBot(i, :), interp1(x, data, atomsBot(i, :)), 'x', 'color', [1,0,0]);
end

% for each TOP layer,
atomsTop = NaN(22, 260);
atomsTopWidth = NaN(22, 260);

for i = 1:22
    i
    data = cutsTop(i, :);
    data = data - smooth(data, 15)';
    [~, idx] = max(data(40:60));
    idx = idx + 40;
    j = 1;
    
    figure(55);
    clf;

    
    while idx < 4096 - 40
        elements  = (idx - 9) : (idx + 7);
        midpt = mean(x(elements));
        y = data(elements)';
        y = y - min(y);
        
        [fitResult, fitgoodness] = fit(x(elements)', y, lorentzFn, 'StartPoint', [max(data(elements)), midpt, 0.03], 'Robust', 'off', 'Algorithm', 'Levenberg-Marquardt', 'TolFun', 1e-10, 'TolX', 1e-10);

        atomsTop(i, j) = fitResult.c;
        atomsTopWidth(i, j) = fitResult.f;

        [~, idxA] = min(abs(x - fitResult.c));
        [~, idx] = max(data(idxA+7 : idxA+25));
        idx = idx + idxA+7;
        j = j+1;
        
    end
    
    figure(21);
    clf;
    plot(x, data);
    hold all;
    plot(atomsTop(i, :), interp1(x, data, atomsTop(i, :)), 'x', 'color', [1,0,0]);
end

%% S2E: Plot the spread of each atomic column. Needs previous section to run first (fits for atomic positions)
figure(3);
clf;

subplot(3, 1, 1);
    [N, edges] = histcounts(atomsTopWidth(1, :), linspace(0.05, 0.2, 100));
    bar(mean([edges(1:end-1); edges(2:end)]), N, 'hist'); 
    xlim([0.08, 0.18]); ylim([0, 20]);
    hold all;
    
subplot(3, 1, 2);
    [N, edges] = histcounts(atomsTopWidth(3, :), linspace(0.05, 0.2, 100));
    bar(mean([edges(1:end-1); edges(2:end)]), N, 'hist'); 
    xlim([0.08, 0.18]); ylim([0, 20]);
    hold all;
    
subplot(3, 1, 3);
    [N, edges] = histcounts(atomsTopWidth(17, :), linspace(0.05, 0.2, 100));
    bar(mean([edges(1:end-1); edges(2:end)]), N, 'hist'); 
    xlim([0.08, 0.18]); ylim([0, 20]);
    hold all;

figure(4);
clf;

subplot(3, 1, 1);
    [N, edges] = histcounts(atomsBotWidth(1, :), linspace(0.05, 0.2, 100));
    bar(mean([edges(1:end-1); edges(2:end)]), N, 'hist'); 
    xlim([0.08, 0.18]); ylim([0, 20]);
    hold all;
    
subplot(3, 1, 2);
    [N, edges] = histcounts(atomsBotWidth(3, :), linspace(0.05, 0.2, 100));
    bar(mean([edges(1:end-1); edges(2:end)]), N, 'hist'); 
    xlim([0.08, 0.18]); ylim([0, 20]);
    hold all;
    
subplot(3, 1, 3);
    [N, edges] = histcounts(atomsBotWidth(17, :), linspace(0.05, 0.2, 100));
    bar(mean([edges(1:end-1); edges(2:end)]), N, 'hist'); 
    xlim([0.08, 0.18]); ylim([0, 20]);
    hold all;

%% S2F: plot spacing between atomic centers. 

figure(7);
clf;
spreadBottom = nanstd(atomsBotWidth(:, 1+10:end-10), 0, 2);
spreadTop = nanstd(atomsTopWidth(:, 1+10:end-10) ,0, 2);
plot((1:22) / 2, spreadBottom * 10, 'o', 'color', [240,72,65] / 255, 'MarkerFaceColor', [240,72,65] / 255)
hold all
plot([0, 12], mean(spreadBottom([2:11,13:22])) * 10 * [1,1], '--');


plot((1:22) / 2, spreadTop * 10, '^', 'color', [77,177,228]/255, 'MarkerFaceColor', [77,177,228]/255);
plot([0, 12], mean(spreadTop(2:end)) * 10 * [1,1], '--');
xlabel('Layer (UC)');
ylabel('Standard deviation of atomic column width (A)');
xlim([0, 12]);
ylim([0, 0.25]);

figure(8);
clf;
spreadBottom = nanmean(atomsBotWidth(:, 1+10:end-10), 2);
spreadTop = nanmean(atomsTopWidth(:, 1+10:end-10), 2);
plot((1:22) / 2, spreadBottom * 10, 'o', 'color', [240,72,65] / 255, 'MarkerFaceColor', [240,72,65] / 255)
hold all
plot([0, 12], mean(spreadBottom(2:end)) * 10 * [1,1], '--');

plot((1:22) / 2, spreadTop * 10, '^', 'color', [77,177,228]/255, 'MarkerFaceColor', [77,177,228]/255);
plot([0, 12], mean(spreadTop(2:end)) * 10 * [1,1], '--');

xlabel('Layer (UC)');
ylabel('Mean atomic column width (A)');
xlim([0, 12]);
ylim([1,1.5]);


mean(spreadBottom([2,3,4,5,6,7,8,9, 11, 12, 13, 14, 15, 16, 17, 18, 20, 21, 22])) * 1000

